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Granular contact force density of states and entropy in a modified Edwards ensemble 

Philip T. MetzgeiQ 

The KSC Applied Physics Laboratory, John F. Kennedy Space Center, NASA 
YA-C3-E, Kennedy Space Center, Florida 32899 
(Dated: February 2, 2008) 

A method has been found to analyze Edwards' granular contact force probability functional 
for a special case. As a result, the granular contact force probability density functions (PDFs) 
are obtained from first principles for this case. The results are in excellent agreement with the 
experimental and simulation data. The derivation assumes Edwards' flat measure — a density of 
states (DOS) that is uniform within the metastable regions of phase space. The enabling assumption, 
supported by physical arguments and empirical evidence, is that correlating information is not 
significantly recursive through loops in the packing. Maximizing a state-counting entropy results 
in a transport equation that can be solved numerically. For the present this has been done using 
the "Mean Structure Approximation," projecting the DOS across all angular coordinates to more 
clearly identify its predominant non-uniformities. These features are: (1) the Grain Factor ^ related 
to grain stability and strong correlation between the contact forces on the same grain, and (2) the 
Structure Factor T related to Newton's third law and strong correlation between neighboring grains. 

PACS numbers: 45.70.Cc, 05.20.Gg, 05.10.Ln, 05.65.+b 



> 

CO 

a 

C 

o 

(N 
> 

m 
O 
m 
O 

o 



-a 
c 

o 
o 



X 



I. INTRODUCTION 

1. Deriving the Contact Force Distribution 

There have been several attempts to derive the gran- 
ular contact force probability density function (PDF) 
for static granular packings, Pp(F), by using analo- 
gies from thermal statistical mechanics 0, H El H 
The interest arises in part because the empirical Pp(F) 
[1 II S IS IH El El El El El El E3 has an ex- 
ponential tail, reminiscent of the energy distributions of 
thermal systems. However, the overall form of Pp(F) is 
not found in thermal systems, generally having a peak or 
plateau near the average force and a non-zero value at 
zero force as illustrated in Fig. (JIJ. 

In contrast to this form, the prototypical distributions 
found in thermal systems are either monotonically de- 
creasing (e.g., the Gibbs energy distribution), or begin 
from zero probability density at the origin before rising 
to a peak (e.g., the Maxwell-Boltzmann distribution). In 
the non-monotonic cases the rising slope is due to the de- 
generacy of energy states. The degeneracies reflect the di- 
mensionality of the system and dominate the form of the 
distribution at weak energies beginning from the origin. 
Since the forces in a granular medium are vector mag- 
nitudes composed from several Cartesian components — 
implying degeneracy in the force magnitudes — this raises 
the question why Pf(F) does not likewise begin from the 
origin Pf(0) = before rising to its peak? Indeed, a re- 
cent model H 0| predicts that it should. The model 
represents a first-principles approximation for key ele- 
ments of the physics and results in a Boltzmann-type 
equation that is solvable. This produces a Pf(F) that 
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FIG. 1: Linear plot of the PDF Pf (/) of the normalized vec- 
tor magnitudes of the granular contact forces resulting from 
Monte Carlo solution of the Mean Structure Transport Equa- 
tion. It has a non-zero probability density for zero force, a 
peak just below / = 1, and an exponential tail with decay 
constant (3 = 1.6. The smooth curve is a fit to Eq. JUJ. The 
log- log inset shows the behavior below f — 1. The dashed 
line is a power law of exponent a = 0.3. These features are 
consistent with experimental and simulation data. 



begins from -Pf(O) = 0, rises to a peak, and then decays 
exponentially. Because of these considerations, the ques- 
tion may be asked whether the empirical observations 
that Pf(0) > is primarily the result of numerical or 
experimental uncertainties: the distribution is in ques- 
tion precisely where the forces are weakest and therefore 
most difficult to model or measure. Perhaps the theory 
provides a clearer view into the fundamental organiza- 
tion of the Density of States (DOS) in this region than 
the empirical methods are presently able to provide. 

It seems to the author that this is not the case for 
two reasons. First, it has been shown that the form in 
the region of weak forces evolves in a predictable way 
as a function of stress and/or fabric anisotropy, which 
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may be induced through shearing [6|. The anisotropy 
dependence probably explains the variations in Pf(F) 
seen among the different empirical studies, in that some 
jammed packings have displayed peaks while others have 
displayed plateaus or monotonic forms. For a packing of 
grains originally in an isotropic state, Pf(F) displays a 
form similar to Fig. ■ As the packing is quasi-statically 
sheared the distribution smoothly evolves to having a 
plateau in the region of weak forces and then on to be- 
coming a monotonically decreasing function with only 
an abrupt change of slope where the peak had previously 
been. After the packing achieves peak shear strength, 
continued shearing reduces the stress anisotropy and 
causes the distribution to retrace its evolution most of 
the way, ending with a small peak again. This behavior 
affects the distribution well above the region of numeri- 
cal uncertainty and cannot be dismissed as the result of 
dynamical or transient forces since the shearing is quasi- 
static. It is difficult to see how this smooth variation of 
forms — including plateaus and monotonic forms — could 
be explained if the finite -Pp(O) > were not real. 

Second, the unique features of the PDF have been ob- 
tained using a wide range of empirical techniques, and it 
does not seem reasonable that all of them are incorrect 
in the region of weak forces. These techniques include 
experiments with frictional grains B, B , numerical sim- 
ulations with friction al g rains 0, Ijjllllll2t or purely 
frictionless grains , an d adaptive network mod- 

els . The simulation techniques have included contact 
dynamics (CD), discrete element modeling (DEM), and 
molecular dynamics (MD) quenched beneath the glass 
transition [l(J, all of which are well-established tech- 
niques. The contact laws in these simulations have in- 
cluded Hertzian, Hookean, and Lennard Jones potentials. 
Simulations have been done with and without gravity 
and under a wide variety of conditions. The transitions 
between the boundary and bulk have been thoroughly 
examined |l5|. The numerical techniques have demon- 
strated the ability to distinguish between distributions 
that begin at the origin and those that do not 1151 . Al- 
though experiments with frictionless emulsions |18| and 
some numerical simulations ^3 have been fitted to forms 
that begin with -Pf(O) = 0, arguably those data would 
be fit as well or better by forms with nonzero Pf(0)- 

The universality of these observations shows that the 
PDF's unique features are not associated with a specific 
type of model or the (non)existence of friction, but are 
fundamental characteristics of static granular packings. 
Because of this, the present paper will proceed with the 
assumption that these observations are correct but have 
yet to be explained. Perhaps the explanation lies in a 
unique generalization of statistical mechanics. Just as 
the DOS for ideal Bose and Fermi gases are organized 
differently than the classical dilute gas and therefore pro- 
duce their own unique energy distributions, so the DOS 
of granular packings may be organized in some unique 
way to produce this distinctive PDF. 

Such a generalization has been taking shape [20I l2ll 



|22T ]. beginning with Edwards' hypothesis |2^| that all 
metastable packings are equally probable in the statis- 
tical ensemble. Another line of progress is based on the 
concept of directed force chain networks [24[ , while others 
aim to understand the distribution of forces beneath a lo- 
calized perturbation or more generally the stress response 
function pB| . an d the phenomena related to jamming and 
unjamming |26L |27|. This paper focuses more narrowly 
upon those models or hypotheses which predict a PDF 
by making assumptions about the DOS in the ensemble, 
including those models which take a random walk in a 
phase space (e.g., the q model) or a PDF space (i.e., the 
Boltzmann transport equation variety), and those which 
directly assume the form of an entropy or other thermo- 
dynamic functional. 

The q model |2^,|2^| may be considered a random walk 
because the set of forces in a single layer of the lattice de- 
scribe a locus in phase space while the random redistribu- 
tion of those forces from one layer to the next (controlled 
by the stochastic q variables) represents a random walk 
through that space. Eventually the walk wanders into re- 
gions of the space having the most probable distribution 
of coordinates. Bouchaud has shown that the sufficient 
requirement to obtain the exponential tail in the q model 
is simply that some grains transmit all their load from one 
hemisphere into just one contact on the other hemisphere 
|22| . This introduced a new way to think about granular 
media: the statistical relaxation of the force distribution 
does not occur dynamically through the time dimension 
as it does in thermal systems; rather it is a necessary 
feature of the internal, layer-by-layer static equilibrium 
relationships, where the spatial dimensions play a role 
analogous to the time dimension for the corresponding 
set of Cartesian components of force [§3 • Several gener- 
alizations of the q model and other lattice-based models 
have been developed [3(]]. Some of these are similarly 
random walks in a non-dynamic phase space, but others 
include explicitly dynamic features to recursively achieve 
organization in the percolating force network. 

In this context it is probably helpful to mention again 
|3l| that the distribution predicted by the q model (2^ 
was not Pp(F), but rather P w (w) where w is the total 
vertical load supported by the grain. Distributions of w 
and F have been occasionally confused with one another, 
especially since w and F become identical in the special 
case at the fiat sides of a container. This has contributed 
to the confusion over the form of Pp (F) . The q model can 
also produce distributions Px{F x ) of the vertical Carte- 
sian components of the contact forces, F Xl but it cannot 
directly predict the vector magnitudes F of those same 
contact forces. The Px{F x ) predicted by the q model 
is always monotonically decreasing, in agreement with 
numerical simulation data |15| . 

Another theoretical model that makes direct state- 
ments about the contact force DOS is the Boltzmann- 
type equation presented by Edwards and Grinev 0, 
mentioned above. In the discussion section, this paper 
shall attempt to reconcile the model with the empirical 
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data. 

Other models include several entropy maximization or 
functional minimization concepts. These methods pro- 
duce elements of the empirically-observed PDFs, but not 
all of their features. The concept proposed by Bagi E10 
deals, like the q model, with Cartesian components. It 
produces the same Canonical distribution as the uniform 
q model. The concept proposed by Kruyt and Rothcn- 
burg 0| deals with contact force magnitudes and predicts 
■Pf(O) = 0, a peak, and an exponential tail. The concept 
proposed by Ngan [j| produces Pp(0) > 0, a peak, and 
a nearly-Gaussian, compressed-exponential tail. Unlike 
Edwards and Grinev's model, these last three are not de- 
rived from first-principles but are hypotheses drawn by 
analogy with other entropic systems. Despite any short- 
comings, all these models provide important insights into 
the nature of the PDF problem. 



2. Organization of the Paper 

This paper is organized as follows. Sec. II will present 
a first-principles analysis of the DOS in a modified ver- 
sion of the Edwards ensemble |2^| . The dynamical be- 
haviors of granular media will be completely avoided so 
that Edwards' hypothesis alone shapes the DOS. It will 
be shown that the DOS is highly self-organized and very 
sparse. Its form depends upon the form of Pf{F) and 
vice versa so that recursion is necessary to solve for ei- 
ther. Maximizing a state counting entropy in this phase 
space will produce the recursion equation which is anal- 
ogous to the Boltzmann transport equation. To eluci- 
date the organization of the DOS, it will be projected 
in the Mean Structure Approximation across all angular 
coordinates before solving numerically. Sec. Ill presents 
the numerical solution of the Mean Structure Transport 
Equation. The results demonstrate the success of Ed- 
wards' hypothesis in that it predicts a form for the PDFs 
in qualitative and quantitative agreement with the exper- 
imental and simulation data, having -Pf(O) > 0, a peak, 
and an exponential tail with a decay constant matching 
empirical observations. Sec. IV discusses the validity of 
the approach and insights into the physics that produce 
the features of Pf(F). Sec. V summarizes the paper and 
points to several unanswered problems and generaliza- 
tions that are needed. 



II. MODIFIED EDWARDS ENSEMBLE 
ANALYSIS 

A. Description of the Particular Ensemble 

Following Edwards and coworkers, this analysis fo- 
cuses upon 2D, amorphous packings of cohesionless, rigid 
grains having the fixed coordination number that makes 
the packing isostatic |2y. The problem shall be further 
idealized, however, by using only smooth, round grains. 



Also, this paper focuses on the frictionless case wherein 
the 2D isostatic coordination number is Z = 4. A method 
has been found to solve Edwards' probability functional 
for this special case. Although this ensemble is highly 
idealized, it is a good starting point because 2D pack- 
ings of cohesionless, round grains that are perfectly rigid 
[II 13 and/or frictionless [IS El E El E ^ are 
known to have force distributions with the same qualita- 
tive features as the more generalized packings and hence 
must be subject to the same organizational constraints in 
the statistics. Therefore they are sufficient to elucidate 
the origin of those constraints in the physics. 

The use of exactly four contacts per grain, however, is 
more idealized than has been achieved in typical numer- 
ical simulations. Nevertheless, it is acceptable because 
that is the average coordination number for 2D friction- 
less packings of round, rigid grains which are isostatic 
[22l l32| , and it will be shown herein that the same qual- 
itative and quantitative features of Pp(F) arise as they 
do in the more realistic simulations. 

Two defining issues for the ensemble are (1) how to 
specify the fabric, and (2) how to apply stress to it. 
Since this paper is concerned primarily with the deriva- 
tion of Pp (F) , and since its form is known to evolve with 
stress and fabric anisotropy under shearing, the ensem- 
ble will sufficiently general to accommodate anisotropy 
in each. On the other hand, this paper does not ad- 
dress the more ambitious problem of stress propagation. 
Therefore the analysis shall not accommodate large-scale 
stress and fabric inhomogeneities that persist in the en- 
semble average. An example of an inhomogeneous case is 
the conical sandpile formed by central pouring, which has 
directional fabric and stress propagation away from the 
center of the pile [3^|. In this paper, only actions such as 
shaking, shearing and compressing are assumed to have 
occurred in the construction history because these pro- 
duce homogeneous stress and fabric states with only sta- 
tistical fluctuations, vanishing in the ensemble average. 

Specifically, stress will be specified as the tensor which 
is a volume average over the entire packing. The source 
of stress will be mechanical at the boundaries to allow for 
the full range of possible states, something which gravity 
alone cannot do. Gravity will be eliminated both because 
it is not necessary and because it breaks a symmetry of 
the ensemble and may thus tend to obscure the organi- 
zational features of the physics. 

This leaves the question how to specify the fabric. Ed- 
wards and coworkers have developed a conjugate field 
theory with the goal of explaining the propagation of 
stresses in granular materials correlated to the local con- 
tact geometry |2l|. In that theory it has been shown 
that two fabric tensors are required to produce the com- 
plete set of stress propogation equations. They have else- 
where developed a thermodynamic theory of compaction, 
in which the relevant specifier is simply the scalar vol- 
ume of the packing [23 • For the sake of simplicity, the 
choice was made there to avoid the full anisotropic treat- 
ment. For the present purposes, something is needed 
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which is less than the full tensorial treatment but more 
than scalar compaction. We will therefore use the joint 
probability density function (JPDF) P4g(8 1: . . . , 84) stud- 
ied in Ref. j34|. This function correlates the contact an- 
gles that share the same grain. It can be collapsed to 
Pe(0), 




...,64)6(8-6,3). (1) 



The use of the uncollapsed distribution is deemed nec- 
essary because the physics of fragile media are grain- 
centered and the intra-grain correlations turn out to be 
the heart of the statistical physics, as shall be shown 
here. This JPDF is not sufficient to relate the pack- 
ing state to the specified stress tensor, however, because 
it tells nothing of the number and size of voids created 
by the arrangement of neighboring grain configurations. 
For the present, the voids may be quantified simply by 
assuming a number of grains per unit distance in a cross- 
section of the 2D packing in each orthogonal direction. 
These quantities along with the JPDF will be specified 
in the ensemble rather than predicted. Evolution of the 
internal state of the packing is beyond the present study. 

Finally, for convenience the idea of "quartered fabric" 
is introduced at several points. It is defined such that 
Pa$(8x, ■ ■ ■ ,84 ) is zero everywhere except where the j th 
contact is on the j th quadrant of every grain in the pack- 
ing. For the specific case of "quartered isotropy," col- 
lapsing the quartered fabric by Eq. (JIJ produces Pg (8) = 
1/27T. This mimics true isotropy but the anisotropic 
quartering is apparent in the JPDF. As in the case of 
non-quartered fabric, P40 enforces steric exclusion. To 
achieve quartered isotropy with steric exclusion in a 
Monte Carlo process it is necessary to weight the dis- 
tribution of attempted angles to emphasize the regions 
close to the edges of each quadrant. Otherwise, steric ex- 
clusion would cause notches to appear in Pg(8) near the 
boundary of each quadrant. The use of quartered fabric 
in this analysis is only to provide insight into the expres- 
sions. It is always possible to write and numerically solve 
them for the more general case, and it was found that nu- 
merical solutions were indistinguishable with or without 
quartered fabric. 

B. The Phase Space 

The locus in phase space of a classical dilute 
monatomic gas completely defines its state. We wish 
to define a phase space for granular packings which is 
similarly complete. A 2D frictionless granular packing of 
N round, rigid grains is isostatic and therefore contains 
2N contacts. The phase space therefore requires at least 
AN phase space axes, half of which define the force on 
each contact and half of which define the contact angles, 
{Fk,8k I k = 1, ...,2N}, which is labeled Si and has 
DOS It is possible to define the ordering of the axes 



so that it is understood which four contacts correspond 
to the same grain and therefore which grains contact one 
another. It will not prove necessary to do so explicitly, 
although this ordering is implicitly assumed to exist. 

Newton's third law (N3L) is automatically satisfied in 
Si, since each contact is represented by only one force 
and one angle axis. However, enforcing Newton's Sec- 
ond Law (N2L) will prove simpler if redundant axes are 
created to account for each contact force twice, one time 
with each grain that shares the contact, {F a p, 8 a p \ a = 
1, . . . , N; j3 = 1, . . . , 4}, where a subscripts the grain and 
(3 subscripts the contact on the grain. This space is la- 
beled S2. In this new space it will be necessary to en- 
force N3L. Again, it is possible to define the ordering of 
the axes so that it is understood which contacts are re- 
dundant to one another and therefore which grains are 
contacting neighbors. It will not be necessary to do so 
explicitly, although this ordering is implicitly assumed to 
exist. 

In the thermodynamic limit N — ► 00 this ensemble 
has Edwards' flat measure, every metastable state being 
equally probable, The DOS in S2 is, 

p m {F a p,e aP } = (5(fabric P 4 e) 




N I 4 \ 4 

x n 9 ^) (2) 

a=l \fl=l I /9=1 

where O is the Heaviside (unit step) function. 

The six constraints which define the accessible regions 
of phase space are described below. 

1. The JPDF for the fabric is specified by the first delta 
function. Actually, there should be a statement relating 
the continuum P4g(8\, 82,83, 84) with the discretized dis- 
tribution of angles at finite N, [ikimn(8ik, ■ ■ ■ , 84-0), but 
the meaning is nonetheless clear. 

2. (and 3.) The Cartesian loads w x and w y on each 
grain will often be called the "supported loads" or simply 
the "loads" . At each locus in phase space the relationship 
exists between these loads and the Cauchy stress tensor. 
(See for example Ref. For simplicity the sum over 
these Cartesian loads is specified. Hence, 

N N 

^ w xa = W x , ^2 w ya = Wy, (3) 

a— 1 a— 1 

where the loads are defined by, 

«fca = + <?*) A Wya = «? + /2, 

(4) 
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and, using non-quartered fabric, 

The operator L a p multiplies the expression by 1 if tt/2 < 
< 37r/2 meaning the f3 th contact is on the left half 
of the grain, else it multiplies it by 0. Likewise R a p, T a p 
and B a p test for contacts on the right, top and bottom 
side of the grain, respectively. For stable grains, w xa = 
w£a — ^"I 11 *! but the hemispheric distinctions shall be 
useful in the analysis. 

4. N3L is satisfied between every pair of contacting 
grains. 

Fyt = -F eC , (6) 

where grains 7 and e are contacting neighbors through 
their 8 th and ( th contacts, respectively. 

5. N2L for static equilibrium must be satisfied at each 
grain individually, 

E^ = 0Vq - (7) 



6. enforces no tensile contacts anywhere in the pack- 
ing, which restricts the DOS to the first "quadrant" of 
the force axes. 

In addition to these six constraints, two missing con- 
straints should be noted: 

1. The above ensemble does not enforce the shear 
stress but relies on the fact that their ensemble average is 
zero and in the thermodynamic limit the fraction of pack- 
ings in which the shear deviates from zero by more than 
some arbitrarily small amount will vanish. The Cartesian 
axes of the packings are taken to be aligned with the prin- 
ciple stress axes so that the off-diagonal elements of the 
stress tensor should be zero. 

2. A cluster of real grains must just touch one another, 
forming closed loops, but in the above ensemble the ge- 
ometric constraints for grains outside the first coordina- 
tion shell have been intentionally omitted. This First 
Shell Approximation (FSA) asserts that only negligible 
correlative information travels all the way around closed 
loops of grains in the ensemble average. In other words, 
the DOS is adequately characterized for the present pur- 
poses by the two-point (intra-grain) force correlations 
and the resulting correlation of loads in neighboring 
grains. Therefore, the geometric closure of force loops 
can be ignored when deriving the statistics of single-grain 
states. There are important arguments supporting the 
FSA and they will be presented in the discussion section. 



C. Phase Space Operations to Quantify the 
Non-Uniformity 

Although Edwards' flat measure is uniform across the 
regions of accessible phase space, the volume of those 
regions is not uniformly distributed across the coordi- 
nates. The program is to change coordinates in a way 
that eliminates delta function constraints from the right- 
hand side of Eq. (J2J), trading the lack of uniformity in the 
constraints for a lack of flatness in the measure. When 
only extensive, conserved quantities remain in the list of 
constraints, then the method of the most probable distri- 
bution may be used, relying on the method of Lagrange 
multipliers to conserve those quantities. 



1. Newton's Third Law 

In this context the term "grain configuration" refers 
not only to a grain's contact geometry but also to the set 
of forces upon those contacts as defined by the locus in 
phase space. The form of §2 itself does not require neigh- 
boring grains to satisfy N3L, and so the vast majority of 
loci include neighboring grain configurations with physi- 
cally unrealizable forces. This mathematical abstraction 
enables the analysis. 

If we wished to neglect N3L then we could proceed with 
the remainder of the analysis, obtaining the hypothetical 
DOS for regions of this space having stable, cohesionless 
grains, write the state-counting entropy and then maxi- 
mize it subject to the conservation of total loads and fab- 
ric. This would produce the most probable distribution 
for all possible permutations of N stable grain configu- 
rations where the grains are mechanically independent. 
However, it turns out that N3L is not negligible: by 
considering instead all the possible combinations (rather 
than all possible permutations) of N stable, independent 
grain configurations, we note that some of these com- 
binations can be mechanically connected into a greater 
number permutations satisfying N3L than can other com- 
binations. Hence, those combinations are the more en- 
tropic ones, the ones which represent the greater number 
of metastable packings in the phase space. Therefore, to 
find the most entropic combination of stable, indepen- 
dent grain configurations, we will map §2 — *■ §3, a space 
where the axes are the same as in §2 except that they 
are not sequenced to represent a particular permutation 
of the grains. Whereas a locus in §2 represents a single 
state (a single packing permutation of a set of grain con- 
figurations) , a locus in S3 represents a set of mechanically 
disconnected grain configurations that may or may not 
be permutable into some number of stable states. We 
shall call the latter an "assembly space" to distinguish it 
from a phase space that identifies every grain's location 
in the packing. The fraction of permutations that satisfy 
N3L will be quantified in this mapping process. 

To verify that a particular permutation of a given 
combination of stable grains {F a p,9 a p} satisfies N3L, 



it is necessary to check every contact in the permuta- 
tion. All permutations of this combination have the 
same JPDF of forces and contact angles, Pfb(F,9) = 
Pfb{F,9 I {F a {j, Q a p}), whatever its form may be. Ran- 
domly choosing one contact from the set of these permu- 
tations, a contact force F a p at angle 9 a /3 therefore has the 
probability PFe(F a p,9 a f3)dFd9 that it will satisfy N3L 
with its neighbor. (The two differentials reflect the fact 
that N3L reduces the solution space by two dimensions 
per contact, thereby taking out the extra dimensions in- 
troduced in §! — ► §2.) The probability that an entire 
grain configuration drawn from this set of permutations 
will satisfy N3L with its four neighbors will be called 

T (-fali ■ ■ • ) Pa4: @al: • • • : @a4 

) d 4 F d 4 6>, written for com- 
pactness as T 2 (F a p, 9 a p) d 4 F d 4 9. It may be written as 
a functional of Pf6, 



to the assembly space, §> 2 — ¥ S3 



T 2 (i 7 Q /3, a [}) — Y\_ ^ Fe (F a 0i Oaf: 



(8) 



{F af 3, 9 a p} = <5(fabric P ie ) $ Wxa ~ Wj ^j 
8 fe w V * ~ Wy] J] T{F a0 , 9 a0 ) d 2N F d 2N 9 



a=l 

I (3=1 



(10) 



The tilde on p indicates that this density is in an assembly 
space. 



This expression treats the contacts on the neighboring 
grains as if they are uncorrelated because this is a pack- 
ing that was drawn randomly from the set of all possi- 
ble permutations, including the ones which are physically 
unrealizable. Therefore there are no a priori correlations 
between neighboring grains; such correlations arise a pos- 
teriori by selecting the subset of packings that satisfy 
N3L. 

Because of this statistical independence, the fraction of 
packings that satisfy N3L for all of its grains is likewise 
simply the product over the probabilities that each of 
the individual grains will satisfy N3L with its own local 
neighbors. (The FSA appears implicitly in this state- 
ment.) However, the product of T 2 (F Q/3 , 9 a0 ) over all 
a accounts for every contact in the packing twice, once 
with each grain sharing the contact. For the cases where 
N3L is in fact satisfied, the double accounting of con- 
tacts will appear as pairs of Pfq factors having identical 
arguments. Hence, the probability that the entire pack- 
ing will satisfy N3L is the square root of that product — 
explaining the use of the square exponent in Eq. |JSJ. The 
fraction of permutations that satisfy N3L is, 



A' 



<fN3L{F a0 ,9 af3 } = [] T{F a ^9 a P) d 2N F d 



2 Nr. 



(9) 



a=l 



This calculation does not handle the boundaries of the 
packing (unless they are periodic), but we are concerned 
with the statistics in the bulk in the thermodynamic limit 
where the boundaries are pushed out toward infinity. 

Now the DOS may be mapped from the phase space 



2. Newton's Second Law 



To quantify the effects of N2L, note that Eq. (JSJ) can be 
used as a many-to-one mapping from S3 — > §4, which will 
have coordinates {w^ a , 9 a p \ £ = x,y;a ~ 1, . . . , N; j3 — 
1,...,4} and is another assembly space, representing 
combinations of mechanically-independent grain config- 
urations. Thus, the mapping reduces the dimensionality 
of the space by two per grain, just as N2L reduces the 
dimensionality of the solution space by two per grain. 
However, the reverse mapping is one-to-one because the 
localized isostacy of the grains determines the four con- 
tact forces when the supported loads and four contact 
angles are specified. Thus, of all the points in §3 that 
map to the same point in §4, at most only one represents 
a stable packing and is occupied, the one which is speci- 
fied by solving Eq. for F aj3 with w x — w l ° {t = w" ght 



and 



= in to P = 



,,bott. 



This system of equations is 



nonsingular except for some precise alignments of con- 
tacts on a grain which we can ignore. The Jacobian of 
transformation for §3 — > §4 is simple to write and is a 
functional of the fabric. Instead, T is re-defined to pro- 
duce the Jacobian with delta functions. 



r 2 (F af3 ,9 a(3 ) 



d 4 Fj[P Fe {F a0 , 



x S (w xa — F a i cos 9 a \ — F a 4 cos 9 a 4) 
x 6 (w xa + F a2 cos 9 a2 + F a3 cos 9 a3 ) 
x 6 (w ya - F al sin 9 al - F a2 sin 9 a2 ) 
x 5 (w ya + F a3 sin 9 a3 + F ai sin 9 ai ) . 

(11) 
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(Note that for simplicity this has been written with quar- 
tered fabric.) With this, the DOS in §4 may be written, 



P {i) {w^ a ,9 a p} = <5(fabric P, 



id 1 



N 

x J]_ T (w xa , w y a , Q/3 ) d 2 N 9 

a=l 

4 

0=1 

3. Cohesionless Grains 



w ya - W y 



(12) 



The product over Helmholtz functions may simply 
be omitted if it is remembered that the Ppg(F,9) fac- 
tors contained within T must be zero for negative ar- 
guments F. On the other hand, it will prove con- 
venient to define ^(F a ^, 9 a p) = Ylp—i 0(F Q( a) where 
Faff = F a p(w xa , W ya ,6 a p). 



D. State-Counting Entropy and Its Maximum 

Randomly drawing packings from the regions of §4 that 
have a specified P2w4e(w x , w y , 9\, . . . , #4) (and hence a 
specified fabric and a specified Ppe(F,8)), the fraction 
of packings in which all grains will satisfy N2L without 
cohesion and satisfy N3L with its neighbors is, 

^{w (a ,9 a/3 } [Pfb] = 

= IL M [Pfb] ■ V(F af} , 9 af} ) d 2N ° 



(13) 



= n» - * - n n ^{w xl ,w V j,9 lk ,...,9 An )[p Fe \ 



x^(w xi ,w y j,9i k , 



d 2N 9 



where Vij k i mn {w xi , w yj , 6 lk , . . . , 4n ) is the discretized 
version of the distribution P2wie(w x ,w y , 61, . . . , 64,), nor- 
malized such that J^i ' ' ' v ijkimn = N. It was ob- 
tained by discretizing the (w x , w y , 6\, . . . , 64) space into 
bins of volume (Aw x ■ Aw y ■ A6q • • • A0 4 ) = (Aw) 2 (A0) 4 . 
Note that in Eq. H13|) , the product in the first line is over 
the grains whereas the products in the second line are 
over the discretized intervals of each of the variables w x , 
w y , 9\, . . . , 64. For compactness we will write, 



<f>{w ia ,9 af3 } [Pi 



F0\ 



n-n[- 



d 2N 9. 



(14) 

To find the most probable P2w4e{w x , w y , 6\, . . . , 4 ) 
that results from the non- uniform DOS of Eq. iJHJ, we 
likewise discretize §4 into cells of volume (Sw x ■ 5w y ■ 



S9 1 ■ ■■80 4 ) N where {5w) 2 {89f = (Aw) 2 (A0) 4 /5, where 
S is a large integer and S » fi... n V (i, . . . ,n). The num- 
ber of cells in §4 which map to a particular set {Vi... n } 
can be estimated by explicit counting, 



u{vi... n ) = II IT 



{S-\ + Vi ... n )\ 
(5-1)! 



E-E> 



(15) 



and in the limit as S — > 00 the estimate becomes ex- 
act. However, because S4 is an assembly space, the axes 
can be relabeled Nl different ways to represent the same 
combination of grains. Removing this physically mean- 
ingless repetition, we omit the factorial of the sums in 
the second line of Eq. lfT5|) , 



u{vi...n) = 11" 'II 



(S-l+I/i,„ n )! 

(s-iy. K.. n )i 



(16) 



where the tilde on u> indicates that this is the "correct 
Boltzmann counting" for an assembly space, in which the 
grains are indistinguishable. 

The number of states fl in the ensemble mapping to 
the distribution {vi... n } is therefore £j{vi... n } times the 
DOS in those cells of §4, 



fiK..»> = n n 



(5-1 



(5-1)! (ui... n )\ 
x AM [T 1 ... n * i ... n ]"*-- d 2N 



(17) 



where we have used the notation of Eq. (|14(l to ex- 
press the DOS. Taking the logarithm, it may be max- 
imized with respect to v p ... u - If we discretize the 
JPDF of fabric P 4 e(0i, • ■ - M -> fi k i mn (e k , . . . , 6 n ) 
such that J2k ' ' ' t^klmn = N, then each value of 
H-kimn^k, I, m, n is a conserved quantity according to the 
definition of the ensemble in which fabric is specified. 
The conservation of W x , W y and ix k i mn are enforced via 
Lagrange multipliers X x , X y and "fkimn respectively. 



i) 



y p...u 



bxQ{vi... n } - X x I E" 'E^ - " Wa;i 

\ i n / 

-\ ^E"'E"' ""'"^ 
E"-E EE"' - 



= 

(18) 



The calculus is performed using Stirling's approxima- 
tion and an expansion of the logarithm in a Taylor series 
of Vi...n where necessary. Taking the limit for S — > 00 
while conserving N and then taking the continuum limit 
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obtains 

P2w4e(w x ,Wy,9 l3 ) = Y (w x ,w y ,6p) * (w x , w y , 6p) 

xG(9p) e- x °> Wx - x * w * (19) 

where the Fabric Partition Factor G(6r) = G(9i, . . . , 9 4 ) 
derives from exp(— 7 P ...«) in the continuum limit. 

Note that and X y should not be confused with the 
decay constants of the exponential tails in the empirical 
PDFs. Most (if not all) of the exponential behavior in 
Eq. I|19l) is contained in the form of T because it is a 
functional of Pfb (F, 9) . It will be shown in a numeri- 
cal solution of the isotropic case using an approximation 
method (later in this paper) that the value of X x = A y is 
approximately zero. This should not be the case in gen- 
eral, however because these two parameters provide the 
only information about stress anisotropy in the equation. 

The Fabric Partition Factor G, along with T and 'i', 
determines the partition of fabric between the (w x ,w y ) 
"modes" . Integrating Eq. (|19f) over w x and w y 

P 4 e(9p) = G{9 fi ) J J d 2 w rfe- x * w *- x v w y 

= G{9 P ) H{9 P ). (20) 
Assuming the standard result of statistical mechanics, 



This can be simplified by taking advantage of symmetries 
in the ensemble. For example, in the case in which fab- 
ric is not quartered so that every contact (/3 = 1, ... ,4) 
is statistically similar, then the summation may be re- 
moved. 

The dependency of T upon the form of Ppg(F, 9) = 
Pf8(F, 9 I 9 a p}) reveals that the DOS in a granu- 

lar ensemble is self-organized (accomplished by the pack- 
ing in its dynamic state as it sought a stable locus in 
phase space) and cannot be given a simplistic a priori 
characterization in a way that is analogous to the a priori 
uniform characterization of a thermal DOS. The form of 
Pfo(F, 9) derives from the DOS non- uniformity and vice 
versa. In principle, this recursion is the unique solution 
for this special case, assuming that correlative informa- 
tion is non-recursive in the ensemble average and that 
Edwards' flat measure over all metastable states is cor- 
rect. 



one form of Pfb will be found in the overwhelming ma- 
jority of the occupied phase space, and so in the thermo- 
dynamic limit we may treat fabric as if it is partitioned 
by G with a fixed form in all of phase space. This factor 
is not a function of w x or w y . However, the partition is 
not an equipartition because of the influence of T and 
The former is variable over the range of angle configura- 
tions within each mode, and the second is a truncating 
factor which limits the range of angle configurations dif- 
ferently for each mode. 



E. The Recursive "Transport" Equation 

The JPDF P2WA8 can be collapsed, 

-1 ' p poo r r r p^tt 

xS[F~Fp{w x ,w v ,9 1 ,...,9 4 )} 
xP 2w4S (w x ,w y ,9 1 , . . . ,9 4 ). (21) 

Inserting Eq. (|19|l into Eq. (|21|l results in a recursion of 
P F e, 



(22) 



F. The Mean Structure Approximation 

The recursion equation can be solved using Monte 
Carlo integration. Efforts are underway to obtain the 
numerical solution, which shall be presented in a future 
publication. For the present an approximation will be 
introduced, simplifying the recursion equation while yet 
providing sufficient accuracy to demonstrate the principle 
organizational features of the ensemble. The approxima- 
tion has value in its own right because it will isolate and 
identify those organizational features. 

The approximation is obtained by projecting the DOS 
in §4 — > §5, where the latter is the subspace {w xa , w ya }. 
This projection is performed by integrating the DOS 
across all the angle axes. For a given pair of val- 
ues (w x i, w y i), the evaluation of T(w x i, w y i, 6 a p) ^ 
Y(w x i,Wyi,6 y s) f°r {Opja {9s}~/ m general. That 
is, certain contact angle configurations {9p} will yield 



1 4 f r°° 

w = iE // d 



0=1 



' 2 W e ~ x * w *~ x y w y 



2tv 



d 4 5(9 - 9 P ) 6 [F - Fp{w x ,w y ,9 ll . . . , 9 4 )} 



G(9 1: ...,9 4 ) f[Q[F 1 (w x ,Wy,6 1 ,...,e 4 )] fff[d 4 F' f[ \P FB [F' s , S )]* 



7=1 



(5=1 



xS[w x - w^(F[, 9 U ...,Fl, 9 4 )} 5 [w x - wi clt ( ■ )} S [w y - w«"( ■ )} 6 [w y - ^ bott '( • )] 

I 
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Region where set of 




^ F 2 ,d2 

FIG. 2: Schematic diagram of 8-dimensional space to illus- 
trate the Mean Structure Approximation (MSA). The MSA 
assumes that the probability for a grain to satisfy Newton's 
third law with its neighbors does not vary over any of the con- 
figurations of the grain having fixed Cartesian loads. In fact, 
the exact probability does not vary too widely for the vast 
majority of those grain configurations. Therefore, the MSA 
should produce a distribution of grain configurations that is 
a good representation of the exact ensemble. The individual 
circles represent the region where random grain configura- 
tions, taken to be neighbors for the grain in question, would 
attempt to apply a particular load on each hemisphere of that 
grain. The intersection is the stable region where the MSA 
applies. 

a greater probability that the grain will be consistent 
with their neighbors than will other contact angle config- 
urations. Therefore, information is lost in the projection 
into the §5 subspace. Nevertheless, this loss of infor- 
mation may not be so great that it blurs the principle 
organization of the DOS. Arguments can be advanced 
to show that, over the distribution of all contact angle 
configurations where the grain is stable, 

Y(w x i,w y i,0 a /3) « T(w x i,w y i) (23) 

for most stable grain configurations. Whereas f is a 
truncating factor in the DOS, defining the region where 
individual grains are stable, T is a scaling factor in the 
DOS, indicating how often particular grain configura- 
tions will occur in the ensemble based upon the prob- 
ability that they can satisfy N3L with their neighbors. 
Eq. (|23[1 claims that this scaling is strongly dependent 
upon the values of w x and w y ; but when varying the 
contact angles it does not vary too much over the major- 
ity of that configuration space. This allows the §4 — > §5 
projection to be simplified. 

The approximation shall be called the "Mean Struc- 
ture Approximation," or MSA, and it is illustrated in 
Fig. It is important to distinguish this from the Mean 
Field Approximation (or MFA), which is useful for ther- 
mal systems but not acceptable for granular packings. 
The reason that the MSA may be adequate where the 



MFA fails is because it preserves the exact intra-grain 
correlation of contact forces by N2L and also the approx- 
imate inter-grain correlations of (w xa ,w ya ) by N3L, both 
of which arc lost in the MFA. The validity of the MSA 
shall be evaluated in the discussion section. 

The most probable p2w(w x , w y ) to occur in the §5 sub- 
space can be obtained directly by integrating Eq. lfl"§|l 
over all angles, 

P2 W (w X) w y ) = e- x * w *- x y w y JJJJ*d 4 G(0 P ) 
x\& (w Xl w yi 9p) Y (w x ,w y , Op) 

(24) 

We wish to simplify this in the MSA. With some manip- 
ulation it can be shown that, 

x S(w x - < ight ) 6 (w y - <p) 
X 5(w x - w l x ch ) S (w y - w^ ott -) 

(25) 

where, 

x 5 « ight - Pi cos 0i - P 4 cos 4 ) 
x S (w 1 ^ + F 2 cos 2 + F 3 cos 3 ) 
x 6 (w y op ~ Fx sin 0i - F 2 sin0 2 ) 
x S (w^ ott - + F 3 sin 3 + P 4 sin 4 ) . 

(26) 

This can be interpreted as the JPDF of attempted loads 
and contact angles that the set of all possible packing per- 
mutations (with the specified Pfb) attempts to place on 
any one grain. The star indicates that this is only a con- 
ceptual distribution, not found in stable packings. It can 
be viewed as a mean-field calculation, where the incoming 
forces have been drawn randomly from the entire set of 
grains in the packing permutations. Its domain is there- 
fore not restricted to the set of forces that would make 
a grain stable. Because of this, the pair (w" ght , w y op ) 
should not be too strongly correlated to (u4 eft , w y ° u ') af- 
ter integrating out the angular dependence, 

PU< sU ,™y 0P ' w x tt > w y 0tt ')~ 

(27) 

PL(< sht , < p )^K e£t , ^ ott ')- 

All the angular content of T is in f 4u?4 g (tojt cm ' , Op ), so 
we make the mean structure approximation, 

PUei^M « P4« m ')/16^ 4 (28) 
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for most stable grain configurations, so that by Eq. Q25[l. 



T(W X ,Wy,9j) W KJ'"!, Wy, Wy)} 2 / 4^ (29) 

for most (w x ,w y ,9j). Using Eqs. O, (23, an d (2HI we 

define 



T(w x ,w y ) = P2 W (w x ,w y )/ATT 2 



(30) 



P2 W {w Xl Wy) may be viewed as a mean-field calculation 
of attempted loads (for packing permutations having the 
specified Pfb), which the half-space of the ensemble at- 
tempts to place on the corresponding hemisphere of a 
grain (for each of the two Cartesian loads), and where 
the mean field includes the unstable regions of the phase 
space. However, it must be emphasized that this is not 
a mean-field calculation of loads actually placed on the 
grains, but rather it is the approximate scale measuring 
how often particular modes will satisfy N3L and therefore 
occur in the ensemble. The validity of the approximation 
depends on the relative weakness of T's dependence on 
the contact angles for most stable grain configurations. 
Using the MSA in Eq. (53) , 



P2w(w x , w y ) 



-\ x w x —X v w 



» W * T(w x ,Wy) 



2tt 



d 4 6 G(6 p )^(w x ,w v ,9f3). 



(31) 



Defining, 

$(w x , w y ) 



2tv 



d 4 Gifip) 



4 

Y[e[F 1 (w x ,w y ,e 1 ,...,9 4 )].(32) 



7=1 



we may write, 
P2w(w x ,w y ) 



T(W X ,Wy) V(W X ,Wy) 

(33) 

Finally, writing the DOS in §5, 
p^{w xa ,w ya } = N\ Y[ Y(w xa ,w m ) ^{w xa ,w ya ) 



xS 



W x 



ME' 



W y . (34) 



We may identify &(w x , w y ) as the "Grain Factor" and 
T(w x ,w y ) as the "Structure Factor." These are the pri- 
mary features of non-uniformity in the DOS. Whereas ^ 
derives from the configuration space of individual grains 
(cohesionless N2L), T derives from the configuration 
space of grains connecting together to form a packing 
structure (N3L). These two factors were so-named be- 
cause their separability (in the MSA) and their roles 



may be considered somewhat analogous to the separa- 
bility and roles of the atomic form factor and structure 
factor of x-ray crystallography. 

The meaning of \& can be illustrated easily through a 
change of variables. We notice that for rigid, cohesionless 
grains there is no inherent force scale and hence stability 
must be independent of the overall scale of the forces. 
Thus it is convenient to change variables, 



w v 



t^xa ~t~ Wy a 



ta 



(35) 



The stability of the a th grain is therefore a function of 
s a and the four contact angles, {0p} a , only. With the 
Jacobian J = t, Eq. H24|l can also be written, 

P st (s,t) = T st (M)* s (s)e- (A * +A * )t/2 - (A *- A « )st/2 (36) 

where the notation has been introduced, 

T st (s,t) = t T[(l + s)t/2, (1 - s)t/2] 

= {w x +w y ) T(w x ,w y ), (37) 



and, 



*.(*) = *[(l + a)t/2,(l-«)t/2] 



(38) 



Note that O in Eq. I132II is insensitive to the scale of Fk 
and cares only whether it is positive or negative, and 
hence the t does not appear as an argument of ^(s). 

In these coordinates Eq. i|32|) may be solved very effi- 
ciently by Monte Carlo integration. This has been per- 
formed for the case of quartered isotropy. In the MSA, T 
does not affect the fabric partition, and hence it is easy 
to find the fabric partition factor. The product of the 
weighting for the quartering bias with the weighting for 
the fabric partition was obtained empirically by adjusting 
as required in a Fourier decomposition to achieve approx- 
imate isotropy Pe(0) « 1/27T. The numerical result for 
that case is fit well by a Gaussian, 



*-(*) 



< 1 



(39) 



with c = 7.9. It is shown in Fig. Q with the fit as the 
dashed curve. This indicates that in the isotropic case the 
volume of a grain's stability space is a Gaussian function 
of the individual grain's load-anisotropy, s. 

In contrast to the simplicity of the above result, the 
form of T depends upon Ppo and hence can only be found 
by solving the transport equation. 



G. The Mean Structure Transport Equation 

Just as Eq. (T5|| can be solved recursively, giving us 
the recursive "transport" equation, so can Eq. Ij33(l be 
solved recursively, giving us the "Mean Structure Trans- 
port Equation" or MSTE. 
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FIG. 3: Grain factor fit to Eq. ffity . 



To develop the MSTE, we convert the load distribution 
of Eq. l(3*3*|l into a contact force distribution. This cannot 
be done simply by collapsing P 2w since it does not con- 
tain sufficient information. However, the variables may 
be changed if we first obtain the joint conditional PDF 
Cfb (F,9 j w x ,w y ), so that, 




FIG. 4: Horizontal cross sections through this plot are the 
conditional PDF for /, the normalized contact force magni- 
tudes ((F) = 0.39 when t = 1). White represents higher 
probability density. The vertical axis represents its depen- 
dence upon the s variable with a fixed t = 1. Varying t only 
rescales /. 



F8 



P Fe (F,6)=JJ dw C F e(F,0 \w x ,w y ) P 2w {w Xl w y ). f C FB , * and T, 

(40) 

This PDF can be obtained easily through the same 
change of variables introduced previously, (w x ,w y ) — > 
(s, i), because Cfb (F,6 | s,t) — t ■ Cfb (tF, 6 \ s, 1) and 
the conditional dependency is reducible to the s variable, 
alone. This may obtained by straightforward integration, 



7=1 



d w 



d 4 G{e p ) 



) 8 F - F 1 (w x ,w 



C F9 (F,e\s,l) = IpJJJJ^d^ G{9 P )5{9 -ft,) 
xj[F-.F 7 (s,lA,---,04) 

4 

ne[F n (a,l,0i,...,04)] (41) 



r,= l 



where only one term of the sum is needed in many 
cases due to the symmetries of the ensemble. This re- 
flects the MSA because it assumes that all grains in the 
same (s,t) mode contribute to the integral according to 
the same weight. It can be found by very easy Monte 
Carlo integration, and the result for the case of isotropy, 
C F e(F,9 \s,l) = C F (F I s, 1)/2tt, is shown in Fig. ©. 

Combining this PDF with Eq. and the definitions 



X JjG^ftWa,,^,^, ...,04,) 

IJ=1 



x e 



c w x — A„u;„ 



P 2 * w {W Xl Wy). 



(42) 



P2 W (w x ,w y ) used in this equation may be obtained a 
number of ways that should be equivalent within the ac- 
curacy of the MSA. Two of these have been used in the 
numerical results and were shown indeed to produce iden- 
tical results to within the statistical precision of the data. 
The first is purely consistent with the MSA, assuming no 
necessity for a priori correlation between the loads and 
the contact angles. Furthermore, it assumes no a priori 
correlation between w x and w y . Correlations arise only 
after throwing out unstable grain configurations. That 
is, it assumes a fixed T over the union of two circles in 
Fig. n °t j us t the intersection of all four (the gray re- 
gion). Imposing then throws out grain configurations 
outside of the gray region. This method is, 



Wu 



jFi cos 6*i + F 2 cos 8 2 , 
F* sin 3 + Fa sin 8a 



(43) 
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(note that all four contacts are treated as if distinctly 
different despite the fact that an x-hemisphere and a y- 
hemisphere overlap in one quadrant), and 

PL(Wx,Wy) = CW^K), (44) 

where 

p p oc p p2n _ 

CK) = // d 2 F // d 2 l[P Fe (Fp,60) 
JJo JJo p =1 

x S(w x — Fi cos 61 — F 2 cos 2 ), 

(45) 

p pOO p p1lT '■' 

^K) = jj d 2 F J J d 2 Y[P Fe (F y ,e y ) 

x S(w y — F 3 sin 6*3 — F 4 sin^ 4 ). 

(46) 

The second method, which will also be used in a Monte 
Carlo solution of the PDFs, attempts greater fidelity to 
the micromechanics by imposing a priori correlation be- 
tween w x , w y and {Bp}. If the MSA is valid, then im- 
posing these correlations should be largely superfluous. 
Comparing the results of these two methods will therefore 
test the MSA in Sec. IV. A. The second method, which 
for simplicity is expressed here for the case of quartered 
fabric, is 

w x = Fx cos 6»i + F 2 cos 6» 2 , 

w y = F 2 sin 2 + F 3 sin 3 (47) 

(note the shared contact F 2 ), and 

p p POO 3 

P2 w3 o(w x ,w yi 6 ) = JJJ^ d 3 F Y[P F e(F y ,e^ 

x S(w x — F\ cos 0i — F 2 cos 9 2 ) 
x 6(w y - F 2 sin 2 - F 3 sin 3 ). 

(48) 

Inserting either of these forms of P 2w into Eq. 142|) 
produces an MSA recursion equation in P F e, which is 
the MSTE. It can be simplified by taking advantage of 
the various symmetries of the ensemble. 

The two different forms of P 2w produce two different 
forms of the MSTE. This is striking because one form of 
P 2w contains {P F e) 3 whereas the other contains (P F e) 4 . 
The ability of these two very different transport equations 
to produce the same P F g depends upon the validity of the 
MSA. 



III. RESULTS 

Here, the following nomenclature is used. The vector 
magnitude of the contact forces are denoted by F, their 



distribution is P F (F), and their mean is (F). The cor- 
responding normalized force magnitudes are / = F / (F) , 
which have a distribution P f (f) = (F) P F (f (F)). The 
Cartesian force components in the x direction are de- 
noted by F x , their distribution is Px(F x ), and their 
mean is (F x ). The corresponding normalized Carte- 
sian forces are f x = F x j (F x ), which have a distribution 

P x (f x ) = (F x )Px(f x (F x ))- 

The MSTE in the previous section was solved in a 
Monte Carlo process for the case of isotropic stress and 
fabric, with one further simplification. It was found that 
X x and Ay were not exactly zero in the MSA, although 
they were very tiny ~ 0.01 so that the exponential fac- 
tors were not exactly unity but were nevertheless neg- 
ligible. Therefore, rather than implementing the expo- 
nential weighting exactly, the forces were simply rescaled 
with a flat factor in each iteration to prevent incremental 
growth. This approach is reasonable because the phase 
space for rigid grains has no inherent force scale, the 
growth was very small, and the growth was balanced in 
the x and y components. Hence, the form of the DOS 
should not be greatly affected by this flat rescaling. 

The MSTE was shown to converge efficiently to the 
same stationary state after beginning from several dif- 
ferent initial distributions. The original work was per- 
formed with Mathematica@ solving for approximately 
5, 500 grains. These results are presented in this paper. 
Ongoing efforts with Fortran demonstrate that converged 
solutions can be found for a million contacts in about 1 
minute on a desktop computer. It is quite easy to ob- 
tain data sets of 10 10 grains or greater, making it possi- 
ble to study joint or conditional distributions of three 
or more variables with smooth statistics using only a 
desktop computer. For some applications this provides a 
tremendous computational advantage over the fully dy- 
namic simulations. 

The Pf (/) resulting from the transport method was 
shown earlier in Fig. (IJ. It has all the key characteris- 
tics of granular contact force PDFs. A fit, shown as the 
smooth curve in Fig. was obtained with the form pro- 
posed for the data from the carbon-paper experimental 
method 0, 

P f {f) = a(l - be~ cf2 y- df . (49) 

Using the values a = 3.28, b — 0.85, c = 1.56, and d = 
1.56, the fit is excellent and is in quantitative agreement 
with the range of values reported from both experiments 
and numerical simulations. It should be noted that here, 
as in most of the empirical distributions 0, El , d is 
suspiciously close to n/2. A plausible reason why this 
value arises under isotropic conditions is provided in the 
discussion section. 

For the special case of true isotropy in which 

P F9 (F,6) = P F {F)P e (e) 

= P F (F)/2n, (50) 

changing variables to Cartesian components F x — F cos 
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FIG. 5: Semi-logarithmic plot of the PDF P x (f x ) of the 
normalized ^-components of the granular contact forces f x = 
F x / {F x ). The smooth curve was obtained from Eq. 1521 . The 
semi logarithmic inset shows the behavior below f x = 1. 



is effected in probability theory by the straightforward 



Px (F x ) = / dfl [°°dF ^^- S ( F * ~ F 
Jo Jo 2vr 



cos 6), (51) 



or by evaluating the inner integral and expressing as nor- 
malized forces, 



Px (fx) 



2 (F) 



d0 Pf (f x sec9) sec (9, 



(52) 



where the symmetries of isotropy were used to reduce the 
range of integration in 9. Numerically integrating this 
[Tlj with the Pf (/) of Eq. (gSJ yields the smooth line in 
Fig. JSJ. It fits the numerical Cartesian component data 
from the transport algorithm (shown in the same figure) 
over the entire range. It has a singularity at f x = and is 
monotonically decreasing as demonstrated in numerical 
simulations [ill [l5T |. It is not purely exponential, the 
two knees being indicative of a summation of n th order 
Modified Bessel Functions of the Second Kind, K n ((3 x f x ), 
functions which result naturally when exponential forms 
are used for Pf (/) in Eq. (52jl. 

The only problem with the fit shown in Fig. Q occurs 
in the region of very small forces, / < 0.2. This is the 
same region in which the form of Eq. H49|) could not be 
experimentally verified due to calibration limits. There- 
fore it is not known whether this is the correct empirical 
form in that region [3(j. A better fit can be obtained 
using another form so that it fits excellently over the en- 
tire range including / << 1. This will be accomplished 
starting with the observation noted above, that the two 
knees in Fig. (JjJJ are indicative of K n ((3 x f x ). These two 
knees appear very dramatically in a rotation of the co- 
ordinates, a rotation which is most easily understood if 
performed manually by lifting the edge of the page to- 
ward the eye and rotating it so that the line of sight is 
parallel to the segments of the graph in Fig. (J5J. The 
fit to Pf(f) will therefore be accomplished by fitting the 
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FIG. 6: (Top) The normalized Cartesian force components f x 
from the Mean Structure Transport Method fitted to Eq. 1531 1. 
which appears to be the natural form. The inset shows the 
behavior below f x — 1. (Bottom) The force magnitudes / 
from the Mean Structure Transport Method fitted to Eq. 1561 . 
The inset shows the behavior below f — 2. These two fits 
analytically transform to one another through Eq. I|52|l and 
Eq. E3. 



natural forms to P x (f x ) and then mathematically invert- 
ing the transformation of Eq. I|52[) . The simplest fit to 
within the statistical accuracy of this data set appears to 
be of the form, 



Px (fx) = Ci anf£K„ (f3xfx) 



(53) 



n=0 



with a = 2, ai = —2, a 2 = 11, and f3 x = 7r/2, and where 
C\ is for normalization. The fit is excellent over the en- 
tire range, displaying all the correct knees and piecewise 
slopes as shown in Fig. ©(top). The shape of the knee 
closest to fx = could be obtained only by including a 
Kq term. This term has infinite probability density for 
fx = 0, but the singularity is very narrow and hence can- 
not usually be seen in a finite set of empirical data that 
has been aggregated into bins of finite width (3?J • 

The transformation integral which is the inverse of 
Eq. (|52|l cannot be deduced by probability theory because 
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f x and 9 are not statistically independent. Therefore, in- 
verting the change of variables to go from (f x , 9) — ► (/, 9) 
is not trivial, even in this isotropic case. Nevertheless, 
the exact relationship can be derived using an approach 
which is equivalent to the mathematics of X-Ray Tomog- 
raphy [^3]. The result is, 

P F (F) = i / d0 P x (F sec 9) esc 2 6, (54) 
F Jo 

or, in normalized forces, 

P f(f) = T^-l [ p * (f sec 0) csc2 6 - ( 55 ) 
\*l J Jo 

This relationship is fascinating because we know that 
F x = Fcos9 and therefore F x < F for all 9; however, 
this relationship computes F in terms of F x — F sec 9 so 
that F x > F for all 9. This says that the probability of 
finding a contact force magnitude F is a weighted sum 
over the probabilities for all the Cartesian components 
F x that are too large to be relevant. Nevertheless it is 
mathematically correct. 

Using Eq. JSJ in Eq. f55|) . we obtain, 

2 

Pf(f) = ^e-efJ2b n (F) n r (56) 

71=0 

with C<x = C\, bo = do, b\ = 7rai/2 + a2, hi = ira.2/2, and 
j3 = p x (F) I (F x ) « (tt/2) 2 . This result fits the numeri- 
cal data from the MSTE excellently over the entire range 
of / as shown in Fig. JBJ (bottom) . It exactly matches the 
finite and nonzero value of Pp (0) = ^C\ao that occurred 
in the numerical data, so we see that the ao term that 
made P x (f x — > 0) infinite is the same bo term that makes 
Pf (0) nonzero and finite. The linear plots of Eqs. I|53|) 
and H56JI are shown in Fig. JJJ) in order to show that the 
curve fits are truly good in the region of weak forces, even 
without the compression of a logarithmic axis. 

Fig. (JHJ shows semi-logarithmically the Cartesian Load 
PDF P w (w) produced by the MSTE, computed for sev- 
eral different rotations of the Cartesian axes. These 
distributions have an exponential tail and a peak near 
w « 1. The near similarity of the rotated plots indi- 
cates approximate rotational symmetry for this nearly 
isotropic model, despite its quartered fabric. The vari- 
ation in the region of weak loads is the result of that 
quartering. In the unrotated axes, wherein the grains 
have exactly two contacts on each hemisphere, we find 
P w (w) — > as w — > 0. We may fit P w (w) in these unro- 
tated axes to an exponential with a power law prefactor, 

P w (w x ) = (^-)' %-/*"-/<«.>. (57) 

If the distribution of F x had been purely exponential and 
if there had been no correlation between adjacent values 
of F x on the same grain, then this should have had values 
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FIG. 7: (Top) Linear plot of the normalized Cartesian force 
components f x from the MSTE fitted to Eq. . (Bottom) 
Linear plot of the force magnitudes / from the MSTE fitted 
□ 



Normalized Load (0 




FIG. 8: Semi logarithmic plot of the Cartesian Load PDF, 
the PDF of the total normalized load borne by each grain in 
the x or y direction unrotated (solid line), rotated 7r/6 radians 
(dotted line), and 7r/4 radians (dashed line). 



of a = 1.0, (3 = 2.0, and (w x ) = 2 (F x ) as in the uniform 
q model. We do find an excellent fit over the entire curve 
using this form, and we do find that (w x ) = 2.0 (F x ), 
but the fit is obtained with the values a — 3 and (3 = 4. 

By comparison, when the Cartesian axes are rotated 
the grains in this model may have 1, 2 or 3 contacts on 
the sampled hemisphere instead of the strict 2 contacts 
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per hemisphere (1 contact per quadrant) that was defined 
for the unrotated axes. The P w (w) for these rotated axes 
are also shown in Fig. (JHJ . They begin with a finite prob- 
ability density for zero force instead of beginning at zero, 
and the finite value is maximized when the rotation is 
7r/4 radians because this is where we obtain the maxi- 
mum fraction of grains having something other than 2 
contacts on the hemisphere. It was found in numerical 
simulations [Til Il5| that when the grains in the bulk are 
segregated into separate populations having one, two, or 
three contacts on one side of the grain, respectively, then 
the Cartesian weight on the grains which support two or 
three others has a P w (w) which does go to zero proba- 
bility for w — * 0. It is the population which supports 
only one contact which has a nonzero P w (w) because the 
load in that case is closely related to Pf(f), which itself 
is nonzero at zero force. Thus, the MSTE results are 
in agreement with this aspect of the simulation data, as 
well. 

The distribution of s and t variables resulting from the 
transport method are fit excellently by 

P st (s,t) = Aco S (^s) (±y e -^e-^ 2 . (58) 

Thus, by comparing Eqs. (|39"|) and with X x = X y = 0, 
the structure factor can be identified, 



T st (s,t) = 




= T s (s) T t (t). (59) 

T t and T s resulting from the transport method are shown 
in Fig. © with smooth curves from Eq. (|59(l . 

IV. DISCUSSION 

A. Validity of the Approximations 

The two approximations which enabled this ensemble 
analysis are the First Shell Approximation (FSA) and the 
Mean Structure Approximation (MSA). Ultimately, the 
quantitative validation of these requires a careful compar- 
ison with numerical simulation data for particular states 
of the stress, fabric, and rheological history, and this has 
not yet been performed. Meanwhile, the qualitative va- 
lidity is already evident as discussed below. 

1. Validity of the FSA 

Beside the constraints which defined the ensemble's 
DOS Eq. J2J), another geometric constraint is needed to 
ensure closure of every "loop" of grains in a packing. 
Without this closure, the chains of contacting grains are 
allowed to branch out ever increasingly in all directions 
and overlap into one another's space. Geometrically, 



then, omitting this constraint does not produce a good 
approximation to a packing. However, it may still be an 
excellent approximation as far as the statistics of single- 
grain states are concerned. 

It has been shown 0] that contact forces on the 
same grain arc strongly correlated with one another. 
There is anti-correlation for contacts closer together than 
A9 w 0.47T radians of angular separation, and a positive 
correlation when the angular separation is greater than 
that. The correlation continues to increase as the con- 
tacts are increasingly distant from one another but still 
on the same grain. The correlation dramatically drops 
immediately thereafter when the distance between con- 
tacts becomes greater than one grain diameter. 

The strong intra-grain relationships make sense due 
to the requirements of static equilibrium of the individ- 
ual grains. Contacts on the same quadrant compete for 
a share of the same load and hence are anti-correlated. 
Contacts opposite one another transmit load through one 
another and hence are correlated. Simplistically we could 
expect A9 — ir/2 to be the crossover point of no corre- 
lation as illustrated in Fig. (|10fl . This is approximately 
correct, and the error is probably attributable to the ex- 
istence of three-grain loops, history-dependent frictional 
effects, and so on. 

Likewise, the sudden drop in correlation after one grain 
diameter of separation is also understandable in terms of 
the local mechanics. It is true that neighboring grains 
share a common contact so that contacts on adjacent 
grains are just two sequential two-point correlations away 
from one another. This induces correlations between 
them. However, these inter-grain correlations should be 
primarily the result of the information contained in the 
sequential two-point intra-grain correlations because the 
lack of cohesion makes the grains otherwise (largely) in- 
dependent. Additional constraints are not found in the 
packing until entirely closed loops of grains are consid- 
ered so that the sequential two-point correlations come 
all the way around the loop back to the original grain, 
again. In 2D the typical closed loop consists of four 
grains, each grain being a vertex between a pair of con- 
tacts that form the loop. The four-point correlation con- 
structed as three sequential two-point correlations going 
the long way around a loop would undoubtedly be very 
weak compared to the single two-point correlation going 
the short way around the same loop, since the short way 
is intra-grain. Hence, the extra correlation information 
imposed going the long- way around the loop must be very 
weak compared to the information already present intra- 
grain. It should therefore be an excellent approximation 
to neglect this additional information and consider only 
the intra-grain relationships in defining the DOS. This is 
the essence of the FSA. 

This is not a rigorous argument because we should con- 
sider the sum of information from all the loops in the 
packing that contain the grain in question, and it is con- 
ceivable that the sum of very many weak contributions 
may be strong. However, due to the randomness of the 
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FIG. 10: Contacts that are approximately 7r/2 radians away 
from one another on the same grain are only weakly correlated 
as illustrated by the closed loop of four grains that allows any 
combination of weak and strong force chains to pass through 
it. If the angles were precisely 7r/2, then the four force chains 
in this figure would be completely independent. 



packing, and the large number of amorphous packings 
that may exist in the configuration space, it is expected 
that the contributions from increasingly larger loops of 
grains will be increasingly decoherent and largely can- 
cel one another. Hence, there is good reason to assume 
that only the intra-grain contribution to the correlations 
is significant in agreement with the FSA. 

If correct, the FSA is an important statement of the 
physics because it fundamentally characterizes the DOS 
and provides deep insight into the organization of the 
physics. In contrast to thermal systems, with granular 
packings it would be completely unsatisfactory to use a 
mean-field approximation because this would throw away 
the structure resulting from the strong two-point correla- 
tions (remembering that these have been observed empir- 
ically). However, by including only this next higher level 



in the approximation, that is, only the two-point correla- 
tions (and assuming that higher correlations exist strictly 
as a sequence of two-point correlations) the maximiza- 
tion of a state-counting entropy and the solution of the 
resulting transport equation produces excellent results as 
shown in the previous Sec. III. The two-point correlations 
therefore appear to be the essence of the physics. Further 
work is needed to carefully test this hypothesis. 



2. Validity of the MSA 

The MSA is important because, if correct, it char- 
acterizes the structure factor as being a functional of 
P2 W (w X) w y ) rather than Ppg(F,9), and this offers the 
possibility to decouple the fabric from the force distribu- 
tions in a way that will help the development of a full 
theory of rheology. In the meantime, pending rigorous 
testing of the MSA, the following three considerations 
are presented to help justify it. 

First, the results produced by the MSA appear to be in 
excellent agreement with the numerical simulation data. 
A focused effort is needed to further test the quantitative 
agreement in specific cases of stress and fabric. 

Second, the values of T have been calculated according 
to Eq. © for the data obtained in the MSTE. The condi- 
tional PDF Py(Y I s,t) was calculated for various fixed 
values of s and t and these are presented in Figs. i|ll|) 
and (|12fl for s — and s — 0.6, respectively. For some 
values of s and t, the ratio T^ max '/T is as high as 3 (or 
greater) and Tf( min )/Y is as small as 1/3 (or lower). This 
means that some grain configurations {si,ti, 9ij} will oc- 
cur three times too often or only 1/3 often enough in the 
MSA ensemble compared to the exact Edwards ensem- 
ble. This effect is most pronounced when t is high and s 
is low. However, high values of t are rare to begin with. 
Furthermore, the distribution for each pair of values (s, i) 
is localized with a clear peak and so the majority of grain 



17 




Log 10 [Y] 



FIG. 11: Distribution of values of T(s,t,9j) for fixed value 
s = and several fixed values of t. (Top) From left to right, 
t = 10 (dashed), 9 (solid), 8 (dashed), 7 (solid), 6 (dashed), 5 
(solid), 4 (dashed), 3 (solid), and 2 (dashed). (Bottom) From 
left to right, t = 1/10 (dashed), 1/9 (solid), 1/8 (dashed), 1/7 
(solid), 1/6 (dashed), 1/5 (solid), 1/4 (dashed), 1/3 (solid), 
1/2 (dashed), and 1 (solid). 



FIG. 12: Distribution of values of T(s,t,9j) for fixed value 
s = 0.6 and several fixed values of t. (Top) From left to right, 
t = 8 (dashed), 7 (solid), 6 (dashed), 5 (solid), 4 (dashed), 3 
(solid), and 2 (dashed). (Bottom) From left to right, t = 1/10 
(dashed), 1/9 (solid), 1/8 (dashed), 1/7 (solid), 1/6 (dashed), 
1/5 (solid), 1/4 (dashed), 1/3 (solid), 1/2 (dashed), and 1 
(solid). 



configurations will have a value of T that is relatively not 
very far from T while being distinctly separate from the 
T for other values of (s,t). These latter considerations 
imply that the MSA does characterize the organization 
in the DOS qualitatively, but more effort is needed to 
show whether it is quantitatively sufficient. 

Third, two different sampling schemes were imple- 
mented as presented in Eqs. 14611 and l]48p. The results 
were identical to within the statistical precision of the 
data, as shown in Fig. I|13|) . This shows that the result- 
ing distributions are insensitive to the existence or non- 
existence of correlations between the Cartesian loads and 
the contact angles, and this is the essence of the MSA. 



B. The Form of the Density of States 

The features of a DOS may be described by two com- 
ponents: the shape of the accessible regions of the phase 
space, and the measure that is used within that space. It 




FIG. 13: Comparison of the curves that were fitted to the 
empirical Pf(f) (large plot) and P x (f x ) (inset) that resulted 
from the mean structure transport method using two different 
sampling methods. In each plot the solid line uses sampling as 
Eq. 1461 with quartered fabric, whereas the dashed line used 
sampling as in Eq. 14811 but with non-quartered fabric. The 
results are statistically indistinguishable, lending credence to 
the mean structure approximation. 
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is possible that the phase space is not equally accessed 
by the dynamics of a real packing as it locates and settles 
into one of the static states. Perhaps this is more true 
for the hyperstatic (frictional) case or for other cases less 
idealized than the one considered here. The form of the 
PDFs would then be a reflection of the shape of the mea- 
sure rather than the shape of the space. Nevertheless, the 
use of Edwards' flat measure produced at least the pre- 
dominant features of Pp(F), and so those features are 
attributable to the shape of the space. The surprising 
repeatability of Pp(F) seen experimentally and in simu- 
lations under many conditions and in many non-idealized 
cases is therefore explained by this fact. 

The rise in Pp(F) to a peak is not due to a degen- 
eracy of F states in the same way that thermal sys- 
tems have a distribution shaped by the degeneracy of 
energy states or momentum magnitudes. In other words, 
if a granular contact force F = y/[F% + F£ + F^] had 
three Cartesian components that were statistically inde- 
pendent, then there would be a volume of phase space 
51(F) corresponding to each value of F such that SI — > 
as F — > 0. If that were the case, then the rise in Pf(F) 
would necessarily begin at the origin. However, since 
that is contrary to empirical observations, this sort of 
Cartesian component degeneracy must not be a domi- 
nant feature in the DOS despite the fact that F is a 
vector magnitude. 

An explanation for this begins with the idea that the 
fundamental unit in a granular packing is not a contact 
force, but a grain, and so the physics of allowable grain 
states limits the space (i.e., there must be a grain factor). 
For the particular case considered in this paper, there 
must be six axes in the phase space of single grain states. 
The 2D stress tensor has two independent principle stress 
values, and so at least two of the six axes represent the 
force state. These may be w x and w y (or s and t). The 
other four axes must convey the geometric information, 
so they may be contact angles. If the space were given 
more axes than this then the density of single grain states 
would be constrained onto a (hyper)surface within that 
space, but we want the states to fill the volume so that 
we may examine the behavior of the volume in the limit 
that one contact force F\ — ► 0. 

A fixed value of F\ defines a 5 dimensional region 
within the single grain space. Its 5D volume is, 

n(Fx) = JJ d 2 wn , (F u w x ,w y ) (60) 

where 

fl'(Fi,w x ,w y ) = JJJJ 6 (steric exclusion) 

x 5[Fi - F 1 (w x ,Wy,9j)]Q(F 2 )Q(F 3 )Q(F 4 ,) (61) 

is the volume of a 3D hypersurface. The integrand of 
this is everywhere nonnegative and for any load state 
(w x > 0, w y > 0) there exist some angles {9p} such that 
the integrand is positive. This is because just three con- 
tacts F2, F3, and .F4 can support arbitrarily high loads 



by themselves regardless of the value of F\ . Therefore 

Q'(F u w x , w y )>0 V w x > 0, w v > 0. (62) 

This fact is demonstrated in 2D frictionless numerical 
simulations where it is seen that a large fraction of the 
grains have coordination Z = 3 and yet support com- 
pressive loads in both axes, w x and w y . Because of this, 
it turns out that in Eq. 1|60|) the integrals in w diverge 
and £1 is infinite for all values of F\. The conclusion is 
that stable grains with F± — > are not confined into a 
vanishing region of the phase space. This is in contrast to 
thermal systems where, for example, p = \J\p\ + Py +Pz] 
can be zero if and only if all its statistically independent 
components become zero so that f2(p) — > as p — > 0. 

There are two key distinctives of the granular phase 
space. First, while contact forces are indeed vectors, the 
stability requirement for the individual grains is so con- 
straining that the components of the vectors cannot be 
statistically independent. Therefore, the DOS cannot be 
uniform in a space defined by the Cartesian axes. The 
degeneracy of vector magnitudes does not automatically 
force Pf(0) to zero. Second, even the magnitudes of the 
contact forces sharing the same grain cannot be statis- 
tically independent. Therefore, the DOS cannot be uni- 
form in a space defined by all the force magnitude axes. 
The vanishing volume of the non-tensile quadrant near 
the origin does not automatically force Pp(0) to zero, 
either. Instead forces sharing a grain are correlated in 
some regions of phase space and anti-correlated in other 
regions, depending upon the the {9p} axes. It is the exis- 
tence of anticorrelation that provides the grains no fewer 
degrees of freedom at F\ = than they have at any other 
value of F\. This will be explained further, below. 

This observation about the phase space is the answer 
why Pf(0) > 0. Edwards' flat measure predicts it, in- 
dicating that the vast majority of metastable packings 
contain a finite fraction of grains with one or more con- 
tacts arbitrarily close to zero force. In fact, we know 
this is correct because every time the stress state of a 
packing is perturbed there is a finite probability that a 
measurable fraction of the grains will tip and rearrange. 
If something in the physics had made the region near zero 
force to be a vanishing fraction of the accessible space, 
then a flat measure in the space would have made tipping 
and rearranging prohibitively improbable. 

Since the volume of phase space does not vanish as 
F — > 0, then what causes the slope in Pp{F) in the re- 
gion of weak forces? The answer is that even though f2' 
is nonvanishing as Fi — > 0, it does get somewhat smaller 
in that limit. This is because contacts on opposite hemi- 
spheres of the grain — say, F\ and F3 — are highly corre- 
lated. When < F\ < (F), then F 3 < over a larger 
region of {Op} than it is when F\ > (F). This was proven 
analytically for a special case in Ref. j3^- Note that 
Eg. 1611 assumes isotropy in the integrand. Weighting the 
integrand anisotropically may provide sufficient general- 
ity to produce either rising or falling slopes in Pp(F) in 
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the region of weak forces, and this may explain its evo- 
lution under slow shearing [(|. 

The reason the Simplest Model of Edwards and Grinev 
predicts P F (0) = is because it treats all the in- 
put forces and angle cosines A 2 as statistically indepen- 
dent. This implies a phase space {Fi,A 2 } with many 
more degrees of freedom than a static grain actually pos- 
sesses. Then, the non-negative domains of all the angle 
cosines ensure that every is positively correlated with 
F, where F = A 2 F\ + . . . + )?( Z -X)F{Z-i)- The only way 
that F can be zero is for all [Z — 1) quantities (AfF,) to 
be simultaneously zero, which is vanishingly improbable 
due to their statistical independence. This is in contrast 
to real grains where the neighboring contacts having less 
than 7r/2 radians of angular separation should be anti- 
correlated. That is, one contact can lift the load off of 
its neighboring contact so that if one contact bears more 
load then the neighbor must bear less. This anticorrela- 
tion allows the grain to be stable with F = while the 
other contacts have nonzero forces on them. That is, the 
grain finds more ways to be stable with zero force on one 
of its contacts than simply by having zero force on all 
of its other contacts. Without addressing the statistical 
independence of the inputs, the model could therefore be 
improved (at the loss of solvability, perhaps) by extend- 
ing the ranges of (A 2 ) to include some portion of (A 2 ) < 0. 
This will account for the range of anti-correlation be- 
tween neighboring contacts. Extending the space this 
way will ensure that Q(F) is nonvanishing for F — ► 
and will result in Pf(0) > 0. This was demonstrated in 
the numerical solution of the MSTE. When the values 
of (A 2 ) were extracted from the numerical data, it was 
found that they had a range —0.5 < A 2 < 1. The lower 
limit reflects steric exclusion, and the upper limit reflects 
maximal separation on opposite sides of a grain. 

V. SUMMARY AND CONCLUSIONS 

The use of the FSA makes it possible to solve the DOS 
based upon Edwards' flat measure in a frictionless granu- 



lar packing of smooth, round, rigid grains with localized 
isostacy. This produces a transport equation that can be 
solved (at least numerically). Solution of this transport 
equation in the MSA was shown to produce the correct 
features for the contact force distributions. 

This success tends to validate Edwards' hypothesis: 
the DOS appears to be dominated by features inher- 
ent to the static phase space, depending solely upon 
the packing's present fabric and the stress tensor. That 
is, the DOS may not be shaped too significantly during 
the physics of the dynamic regime before the packings 
achieve static equilibrium. 

The need for further work is apparent. First, the two 
approximations have not been adequately validated. The 
quantitative results should be compared with simula- 
tions of rigid, frictionless grains with carefully controlled 
stress states and carefully measured fabric. This has not 
yet been performed because most studies have either in- 
cluded gravity or not reported the stress state or fabric. 

Second, solution of the transport equation without the 
MSA is being developed. Those results compared against 
the present study will be an important test of the MSA. 

Third, the analysis should be extended and numerical 
results presented for more general cases. The case with 
anisotropic stresses and fabric should demonstrate the 
qualitative evolution of the PDFs under shearing. This 
work has begun, and the initial results are hopeful. 

Fourth, the forms of the functions that fit the numeri- 
cal data for P s t(s, t) are tantalizingly simple. If the cause 
of this can be identified then a completely analytical so- 
lution to the MSTE may be possible. 
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